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SUMMARY 

In this study the program named FANS-3D (Finite Analytic Numerical Simulation - 3 
rtfmMtcionan is presented FANS-3D was designed to solve problems of incompressible fluid 
StiJbKlXUcr. I. solves problems wish conduction »d l convecnm, 
modes of heat transfer in laminar flow, with provisions for radiation and turbulent can 

solve singular or conjugate modes of heat transfer. It also solves problems in natural comecnon, 
using Ae Boussinesq^ approximation. FANS-3D was designed to solve heat 
inside one two and three dimensional geometries that can be represented by orth g P 
™Cvterian coordinate system. It can solve internal and external flows using appropriate 
boundary conditions such as symmetric, periodic and user specified. 

introduction 

The program FANS-3D solves one, two and three dimensional fluid flow and heat transfer 
problems that involve conduction and convection modes of heat transfer in . . 

Lunar flow, with provisions for radiative heat transfer and buoyant and ^ b ^ t fl ^ S f I ^ S0 
solves problems in natural convection using the Boussinesq approximation. Uangtas feati^ 
Ae proLm may also solve mixed natural and forced convection problons. Furthermore, «can 
Lve indtodual modes of heat transfer as well. The program FANS-3D solves any 
that can be represented by orthogonal planes in a Cartesian coordinate system. can 

solve internal and external flows using appropriate boundary conditions^h^ s>romemc, 
periodic and user specified. The program is designed to have the same performance in all 

directions, in this wiy, any problem can be solved in the most . 

The program FANS-3D is based on the 19-point Finite Analytic Method. It uses the 
SIMPLE? iterative method suggested by Van Doonnal and Raithby (1984) to solve the coupled 
Navier-Stokes equations. The discretization of variables is done following a new stagger gn 
layout. The resulting system of algebraic equations is solved by different methods, including 

^’^e progl^M^Tr^dSes; a graphics and a computational. The graphics module, 
named GRAPH3D, was written in FORTRAN 77 for Apollo workstations. This 
displays the geometry of solution and the results in the three-dimensional space. These are 
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presented in the form of velocity vectors, profiles, and color contours with shading of any 
variable These graphic results can also be sent to monochrome and co or P™ ' . , . 

Uc computational module of the pro^am is divided in two pans^The first, which , » 
accessible to the user should be modified according to the problem. The second part is fixed 
Jdoelnot “Lr intervention. These two parts must be bound together to create the 
eonmutational^un' file Both parts of the computational module are written , n standard 
SEES? 77 W m this form the program can be easily poned to almost any machme. 

THE STRUCTURE OF THE PROGRAM FANS-3D 

As mentioned above the program FANS-3D is based on the 19 Point Finite Analytic 
Method. The basic idea of this method is illustrated here Details of its derivation are given by 
Bravo et. al. ( 1 99 1 ). This method is derive for the general transport equation 

+ + 4>" -D tyl + 2A (J)j: + 2B tyy + 2C <!>; — Sh . .. 

using the Jaiytic fotaioos of simplified fotms of it. To sun this eqmnonts ^'o^ly 
linearized in the three-dimensional element shown in figure 1 . To this effect, the coefficients A 
^^,d the source term S are assumed consunts and equal to their values at the center of the 

element, i.e. lTi 

A_ + <h n , +4>» = Dp$t + 2Apbx + 2Bpfyv + 2Cpb;-S p ..... . V 

Thep subindices mean that these coefficients are considered constant inside the 
equal to their values at the center For example, if + is the u velocity ^mponentin a lunmar 

flow and R the Reynolds number, then A p = Ru^,B -Rv p /2,C =Rwy 2, p 

Equation (21 is then solved in the planes x=0, y=0 and z=0 and shown in figure 2. These 
two-dimensional solutions are combined to obtain the three-imensional 19-point finite analytic 

scheme. Details of this process are given by Bravo et. al. (1991 ). ... 

The solution of the coupled Navier-Stokes equations present an additional inconvenience, 

there is no clear equation for pressure. To solve this problem many methods have been 
developed. Notable examples are SIMPLE (Semi-Implicit Method [or Pressi^Lmked 
Equations) of Patankar and Spalding (1972), SIMPLER (Paumkar, 1980), Sl^LEC O'M 
Dermal and Raithby, 1984) and PISO (Issa et. al., IWg. Most 
known as pressure correction methods. The program FANS-3D uses the SIMPLER and 
SIMPLEC methods. The discretized system of equations when usedwnth the Finite Analytic 

Method can be found in the works of Bravo (1987) and Aksoy (1989). 

Another major difficulty in the implementation of the numerical schemes t0 
fluid flow problems is the choice of a proper computational grid. Clearly, it 
if one could discretize the governing equation using a gnd system that P lac * s 
variables scalar and vector, at the same physical location. Unfortunately, the use of such a 
nonstaggered grid system with a primitive variable formulation of * c 
has been shown to produce nonphysical oscillations in the pressure field (Pa * 
remedy to this problem is the use of a staggered grid system; first introduced by Harlow and 
Welch y (1965). This grid distribution was used successfully m many codes anditissttH the most 
prevalent grid arrangement. There are many advantages of this type of staggered gnd 
LrangemSJ. For a typical control volume, this discretized continuity equation contains the 
differences of adjacent velocity components. The discretized gradient of P r «^ e ™ 
momentum equation also contains adjacent pressure values. This 

occurrence of a wavy pressure and velocity fields in the numerical solutions. In the stagg 
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occurrence of a wavy pressure and velocity fields in the numerical ^ force 

grid arrangement, the pressure difference between two adjacent gnds 
for the velocity component at the cell face between these gnd points. Besides, he mass 
rate across each cell can be calculated without any interpolation for the relevant velooty 

arrangement of the grids also eliminated the need of specifying the 

PrC Xw" S^t has its disadvantages too. In this arrangement, 
there^etwo^ distinct cells foMhe application of the finite an^iemefo^ to Ae 
two-dimensional momentum equation as shown in figure . *s trap > 
two sets of finite analytic coefficients. In the three-dimensional problem °* 

coefficients must be evaluated for the momentum equations, one for each vdocity ' 

This increases the CPU cost and memory requirement. The case of staggered gnds also 
increases the difficulties in programming since different geometnes parameters raus 

To solve these problems Rhie and Chou (1983) introduced a non-staggered grid 
arrant S method all variables are evaluated at the same locatior ° fcach 

control volume The pressure gradient in the momentum equation is sull evaluated by 

“ «» "odes. 

eauation new velocity components on the volume faces are evaluated. These velocity 

otoiL by an interpolation scheme based on the mott^tum ^uanon^ Tbts 
mJrt apparently devised by Rhie and Chon- (1983) was M- > . mdSctad 
(1988). They called it the pressure-weighted interpolation method (P^lM). l hey repo 
the PWIM scheme predicts physically unrealistic velocities in regions o rapi y varyi 

m fim mccessM application of a non-staggered grid arrangement to the finite tmalytic 

method was done by Aksoy (1989). Their rae ' h ^' ^ c ^ ^9«) m^*!d, but with 
Interpolation Method), is similar to the original Rhie and <^w (1983) FW1 M im w. 
a different interpolation scheme. The advantages and disadvantages ofMWIM are simiwr o 
PWIM method There is only one set of FA coefficients to be 

equations, reducing memory and computational time. However, using **“*.’ *, is 

velocity components are also obtained in regions of strong pressure gradien . pro 

^ r^obTei^l^dbove, for die staggered and «m-staggered s^ a 
new scheme is used in the program FANS-3D. In this new scheme, a staggered gnd system is 

used, but only one set of coefficients is evaluated- ^le new 

arrangement, but it also uses the main concepts of the MWIM and PWTM ®<*h • 

scheme evaluates only one set of coefficients at the center of each control volume 
Coefficients at the nodes of the velocity components are obtained bylm<^n^on 
of Ae coefficients obtained at the centers of Ae control volumes. . faKS . 3 d 

given by Bravo (1991). This method was folly tested and it is Ae meAod used in F 
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SOME APPLICATIONS OF THE PROGRAM 


Conjugate Heat Transfer in Laminar Flow Between Par allel Conducting Plates (Extended 
Graft Problem) 


The geometry of the problem is represented in figure 5, where fully developed flow 
between two infinite parallel plates is moving from left to right. The domain has a total length 
of 10H with inlet height of H. The parallel plates have a thickness equal to 0.5H and are kept at 
constant temperature T, on their external surfaces. The flow enters at uniform temperature T f 
The conductivities of the solid and fluid are k t and The conductivity ratio R k ** k/kj 
determines the temperature in the solid fluid interface. If the conductivity ratio R t is very high, 
the temperature on this interface is very close to the external surface temperature T w . In this 
case, the problem is reduced to the original Grata problem (Bravo, 1987). In general, the 
distribution of temperatures is determined by the peclet number Pe, and the conductivity ratio 


**■ 

The Gratz problem, when the temperature on the solid fluid interface is kept constant, was 
solved analytically by Prins, Mulder and Schenk (1949) and later numerically by Bravo (1987). 
The solution for the general conjugate heat transfer problem was later obtained by Mori, 
et.al.(1989). To solve the conjugate heat transfer problem, Mori et.al.(1989), represented the 
interfacial temperature distribution by infinite power series. They solved the governing energy 
equations for the solid and fluid domains using this temperature distribution as boundary 
condition. The Nusselt numbers, obtained by this procedure, were also presented in the form of 
infi nite series. The major sources of error in their analysis are die truncation error during the 
evaluation of this infinite series and the simplification of the energy equation in the fluid flow 
domain. They assumed negligible diffusion in the longitudinal direction. 

To compare with the results of Mori et.al. a Peclet number of 500 was selected. Four grid 
sizes of 10x3x9, 20x3x20, 40x3x40 and 80x3x80 along the x, y and z directions respectively 
were used. Insulated boundary conditions were applied on the solid inlet and outlet boundaries. 
The outlet boundary condition for the flow domain was considered fully developed or dT/dx = 
0 . 0 . 


The dimensionless temperature, for this problem, is defined by 

0 Tf- T w . 


(3) 


With this definition the dimensionless external surface temperature 0 W is equal to 0.0 and the 
dimensionless inlet temperature 0/ equal to 1 .0. The local Nusselt number Nu z is defined by 

H q f 


Nu x - 


where 


k/J fo - T m ) 


(4) 


q { ~ heat flux at the solid fluid interface 
= temperature at the solid-fluid interface 
T w = external constant surface temperature 
7}= uniform temperature of the incoming fluid 
T m = mixed mean temperature defined by T m = 
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With q, = -*/§ , and using .he inlet height H as reference length, equation can be reduced to 
its dimensionless form ^ ^ 

Nu ‘ = 6/0 - ar 

wtKre e /m = im 1, Ojt, the dimensionless interface temperature and »• the dimensionless 

“deprogram FANS-3D finds the temperature 6 at each nodal point. The Nusselt number 
is evaluated from this temperature distribution and equanon (5). 

Figure 6 compares the Nusselt number obtained by with the vjdues 

obtainedby Mori et. al. (1989) and discussed above^For £ thepmgram 

an excellent agreement between the resu s ® values riven by Mori et. al. are under the 

FANS-3D, for x values ov « 0.2. e « ow x _• - ^ Nussc i t nU mber at the entrance must 

ones given by the P™?™ 111 FA ^ S ' 3D , , . th CTam xhe lower values given by Mori 


qute wllvdth the corresponding ones 8’''.° R^ie^j^dytic metitc^ 

values computed by Bravo (1989) ustng the .9 pom, Itate r^y« 

(these values are not shown in the figure). The program FANS-3D provioes 
reasonable and agree quite well with previous computations. 

rnnfnpate Transfer in ■ CtHPPfftl H*** EKltM KfiE 

offset-fin heat exchangers (Kays and London, 1984). To me mp ^ The flow 

exchanger is beyond the capacity of any computer and a an^lvfit^ « 

between the finned plates is three-dimensional and very cmnplex °' h "! 

transfer problem is aggravated kwtk weuse.be 

dimensional conjugate heat transfer problem. To ob^ a ^unontom p S w 
concept of fully developed flow extended to these geometries ficldbwome 

1977) In this case, the velocity, a reduced pressure and a reduced 
periodic after some entrance length. The reduced temperature is defined by 

. T\x,y,z)-T w (6) 

0(x,}\z)— T —T . 

where r_ is the bulk tmpora^ure a, any - ™ S "sh^“ 

flow inside ducts of uniform cross section (Chap , ' , . ^ neriodic length L are 

temperature profiles at successive streamwt^ locations ^amt ^ 
assrnn^ similar. Using this idea, it is possible to simplify 4e geomdjy r P An enlarged 

reduce the computational domain to the one *1 front aid back walls are 

view of this small domain is shown in figure 8. The temperature in the 
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assumed constant and equal to Tw. The separation between the parallel plates is considered the 
reference length equal to L. The position of the fin and its dimensions are given in the same 
figure. The flow enters the domain in right hand side and leaves the domain in the left hand 
side, as shown by the arrows. The top and bottom surfaces are planes of symmetry. 

Assuming constant properties, the temperature field and velocity fields are decoupled and 
their solution can be obtained independently. The governing equations for fluid flow are of the 
same type of the standard transport equation (1) and the program FANS-3D can be used without 
any modification. 

To solve for temperature we define a new variable <|> given by 
x/ r ,. _ v 7\x,y,z)-T w 

y> z ) = — T z f w — 

where T bi is the inlet bulk temperature and T w the external wall temperature. Using this 
definition and equation (6), the profiles of the dimensionless temperature at the <j) 0 at the outlet 
and <|>i at the inlet are of similar shape with the outlet condition being a constant times the inlet 
condition, or 


4*0 = (8) 

where is the bulk exit temperature. To determine the profiles of temperature at the inlet and 
outlet, an iterative procedure is required. We start assigning the temperature at the inlet 

= l .0, then the temperature on the outlet <|> 0 is obtained by solving the governing 
energy equation. Using the values of 4» 0 on the outlet the constant <(>/*, is calculated by 
numerical integration at the outlet boundary. Finally, the temperature at the inlet is updated. 
This process is repeated until convergence is attained. All these steps are automatically 
performed by FANS-3D. 

Figure 9 shows profiles of velocity component u along the y and z directions at the inlet 
and outlet planes, as obtained by the program FANS-3D. In this figure we can recognize the 
periodicity in the velocity distribution and the three-dimensional character of the flow. The flow 
not only moves up and down due to the presence of the beam, but also moves to the left and to 
right. 

Figures 10(a) and 10(b) show the temperature profiles at three locations, at the inlet, on the 
beam and at the outlet. In these figures we can observe that the minimum dimensionless 
temperature <J> is located on the lateral walls. The maximum temperatures are on the plane of 
symmetry. The temperature inside the fin, also shown in this figure, is almost constant in each 
cross section, but increases towards the center between die plates. 

To study the effect of the fin, the pressure drop and the rate of heat transfer in this geometry 
can be compared with the ones for a fully developed flow between parallel plates. For example 
the pressure drop in the element of the heat exchanger, as computed by the program FANS-3D, 
is equal to -0.4565. For the equivalent situation, but without the fin the pressure drop is only 
-0.1200. Therefore, the increase in pressure drop due die fin is 280%. 

The total increase in energy content, when the fluid moves from inlet to exit is given in 
Table 1. This table shows the energy increase in the element of heat exchanger for water and air 
and compares it with the one in a fully developed flow between parallel plates. The effect of the 
fin is an increase in the global rate of heat transfer of 68.7% for water and 54.2% for air. An 
additional computation with an infinite conductivity of the solid showed a less than 0.2% change 
in the rate of heat transfer. Thus, to improve the design of the compact heat exchanger, it is 
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fet m stored distribution of fins can be more affective than the regular one. 

CmUMilt Hea t laMfo I" Electronic Modu les 

In this section the combined effect of convection, and radiation heat ttansfa is studied in 
an airayof electronic chips displayed in figure 11. The chips 

Method and the radiation analysis using the discrete ordinates ^*^Wn^iSon and 
di , con) Results showing the effects of convection alone, combined radiation ana 

Son. rL R pmsence or no! of a tadif vely pamc.panng — ^e = ed fo 

To cimnlifv the solution fully developed penodic flow is assumed (see discussion ior 

prevL p^Using this concept of pmodicity 

of the planes of symmetry, the problem can be reduced to the ep \ 2 ?Tte fluid 

1 1(b) and 11(c). The dimensions in centimeters of this element are gi gur 

V£ «^ * *< ««* is A11 

(air) with constant properties. The flow is considered laminar. 

( TheNavier-Stokes equations are decoupled from the energy oquaoon and *e fluid field can 
, „ lvftd indenendentlv The solution was obtained using the program FANS-3D with the 
19-point Finite Analytic method. Once this solution was performed, the energy 
^IvTilingthe same method explained above for the compact heat exchanger The oou^mg 

models arc interdependent through the temperature field and the wall Tfac 

The results for the flow field compuution are shown m figures ^ ^ tfac 

velocity profiles for the u component are shown in figure 13 at two locatio , 

££ between the blocks. From these figure, we can predate 
blocks and the upper plate is similar to the velocity between two plates, or Poisemlle : flo 
Reynolds number of 1 00 the velocities between the blocks are quite small, gure *P 
Sty vS« on the back plane of symmetry. In this last figure, at the ^ ^AeW^s^e 
flow s rotating counterclockwise. Figure 15 show velocity vectors on a plane l^ ^meen 
Sie blocks in the x direction. This figure shows that some fluid is entrained from the top of the 

blocks and is transported to the lower sides. r- . i c 

Profiles of dimensionless temperature <)> on three different planes are shown in 
and 1 7 The temperature of the blocks is zero and the bulk inlet temperature isequ • 

Figure 16 shows these profiles when convection heat transfer alone is considered. P« 
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the profiles are perpendicular to the top and bottom insulated plates. In this figure we can also 
appreciate that the fluid moving between the blocks has a temperature very close to the 
temperature of the blocks. Specifically, the fluid near the bottom wall has almost the same 
temperature of the blocks. The maximum dimensionless temperature (coldest dimensional 
temperature) is on the top insulated wall. 

Figure 17 shows dimensionless temperature profiles when convection and radiation are 
included. The temperature profiles are quite different from the ones discussed in figure 16. The 
maTimum value of the temperature is in die center region between the top plate and the blocks. 
The top and bottom plates have an intermediate temperature between the temperature of the 
blocks and the temperature in the center region. The higher temperature of the top plate and the 
lower temperature on the bottom plate, compared with the convection case alone shown in 
figures 16, is because of the radiation effect. The top plate receives a net radiation coming from 
die blocks and the bottom plate. Because this plate is insulated, this arriving net radiant energy 
is dissipated by convection, which is indicated by a positive slop of the profile at this point. The 
<}> profile has a maximum in the center region and decreases as we move closer the bottom plate. 
However, close to the bottom plate this profile increases again. This increase indicates that heat 
is transfer by convection from the fluid to the bottom wall. Because this wall is insulated, the 
same amount of energy is irradiated to the top wall. 

An energy balance when convection-radiation heat transfer was considered shows an 
increase on the rate of heat transfer of 49.1% compared with the computation considering only 
convection heat transfer. In this problem then radiation is a very important and must be 
considered. 

This final problem of electronic modules is an example of a three-dimensional 
conduction-convection-radiation heat transfer problem. Although conduction was not explicitly 
discussed, the problem was solved assuming infinite (TO 50 ) conductivity in the blocks. In this 
situation the blocks assume a constant temperature everywhere resembling an isothermal body. 
The convection and convection-radiation results show the potential of the program FANS-3D in 
the simulation of complex three-dimensional problems that include all modes of heat transfer. 

CONCLUSIONS 

This study shows the solution of complex three-dimensional problems that included 
conduction convection and radiation modes of heat transfers with the application of the program 
FANS-3D. The study of these problems also shows different types of boundary conditions from 
die simplest boundary when the values of the variable are assigned on the boundary, to 
symmetric and periodic boundaries. The solutions were presented in the form of vectors and 
profiles given by the graphics part of the program FANS-3D. This giaphic program also 
displays contours that include shading. The difficulties in reproducing these colors do not allow 
the inclusion of these pictures in this work. 
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Table 1 Comparison Between Energy Increase In tbeElcment of Heat Exchanger 
and Energy Increase in Flow Between Parallel Plates 



Water 

Air 

Pr 

1.76 

0.704 

k(W/m°C) 1 

0.6775 

31.27x10*3 

AEwith fu/Cl 

0.07075 

0.15658 

AEnarallel olatcf^ 1 

0.04194 

0.10157 

Inc. % 

68.7 

54.2 



Figure 1 Finite Analytic Element 



Figure 2 19 -point FA Method 


195 




Figure 3 FA Cells in a Staggered Arrangement 



Figure 4 Evaluation of FA Coefficients 
in the New Method 
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Figure 5 Domain of Solution for the Extended Giatz Problem 



Figure 6 Local Nusselt Number at Different Locations 



Figure 7 


(c) 


Heat Transfer Analysis in a Compact Heat Exchanger 
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(b) 


Figure 10 Temperature Profiles on Three Different Planes 
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Figure 12 Domain of Solution of the Electronic 
Components Shown in Figure 1 1 



Figure 13 Profiles of Velocity Components U 
3-D View 
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Figure 14 Velocity Vectors on Plane of 
Symmetry y = 0 



Figure 1 5 Velocity Vectors on Plane x - 1.5 
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Figure 16 Dimensionless Temperature Profiles 
Conduction-Convection Problem 



Figure 17 Dimensionless Temperature Profiles 
Conduction-Convection-Radiation Problem 
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